#
# BAM_CHESUS2017.r
#
#
library(foreign)
library(plyr)
library(dplyr)
library(questionr)
library(statar)
library(ggplot2)
library(viridis)
library(rjags)
library(apsrtable)
library(mokken)
library(lattice)
library(GGally)
library(reshape2)
library(MASS)
library(gridExtra)
#
#
CHESUS <- read.dta("CHES-US_2017.dta", convert.factors = FALSE)
#
# Table A1
print(aggregate(CHESUS[,4:10], list(CHESUS$country), function(x){round(mean(x, na.rm=TRUE),2)}))
#
# ALL PLACEMENTS
#
US.expert.placements <- t(CHESUS[CHESUS$country=="US", 5:46])
EU.expert.placements <- t(CHESUS[CHESUS$country!="US", 5:46])
#
nstimuli <- nrow(US.expert.placements)
df.plot <- matrix(NA, nrow=nstimuli, ncol=4)
for (i in 1:nstimuli){
df.plot[i,1] <-t.test(US.expert.placements[i,], EU.expert.placements[i,])$estimate[1]
df.plot[i,2] <-t.test(US.expert.placements[i,], EU.expert.placements[i,])$estimate[2]
df.plot[i,3] <- t.test(US.expert.placements[i,], EU.expert.placements[i,])$p.value
}
#
df.plot <- data.frame(df.plot)
rownames(df.plot) <- rownames(US.expert.placements)
names(df.plot) <- c("meanUS","meanTA","pvalue","order")
#
df.plot <- df.plot[order(df.plot$pvalue),]
df.plot$order <- 1:nstimuli
#
rownames(df.plot)[rownames(df.plot)=="partyaecon"] <- "Party A (econ)"
rownames(df.plot)[rownames(df.plot)=="partybecon"] <- "Party B (econ)"
rownames(df.plot)[rownames(df.plot)=="partycecon"] <- "Party C (econ)"
rownames(df.plot)[rownames(df.plot)=="partyasocial"] <- "Party A (social)"
rownames(df.plot)[rownames(df.plot)=="partybsocial"] <- "Party B (social)"
rownames(df.plot)[rownames(df.plot)=="partycsocial"] <- "Party C (social)"
rownames(df.plot)[rownames(df.plot)=="partyasocial"] <- "Party A (social)"
rownames(df.plot)[rownames(df.plot)=="partybsocial"] <- "Party B (social)"
rownames(df.plot)[rownames(df.plot)=="leftrightdem"] <- "DParty (leftright)"
rownames(df.plot)[rownames(df.plot)=="leftrightrep"] <- "RParty (leftright)"
rownames(df.plot)[rownames(df.plot)=="leftrightbernie"] <- "Sanders (leftright)"
rownames(df.plot)[rownames(df.plot)=="leftrighttrump"] <- "Trump (leftright)"
rownames(df.plot)[rownames(df.plot)=="econdem"] <- "DParty (econ)"
rownames(df.plot)[rownames(df.plot)=="econrep"] <- "RParty (econ)"
rownames(df.plot)[rownames(df.plot)=="econbernie"] <- "Sanders (econ)"
rownames(df.plot)[rownames(df.plot)=="econtrump"] <- "Trump (econ)"
rownames(df.plot)[rownames(df.plot)=="socialdem"] <- "DParty (social)"
rownames(df.plot)[rownames(df.plot)=="socialrep"] <- "RParty (social)"
rownames(df.plot)[rownames(df.plot)=="socialbernie"] <- "Sanders (social)"
rownames(df.plot)[rownames(df.plot)=="socialtrump"] <- "Trump (social)"
rownames(df.plot)[rownames(df.plot)=="immdem"] <- "DParty (immig)"
rownames(df.plot)[rownames(df.plot)=="immrep"] <- "RParty (immig)"
rownames(df.plot)[rownames(df.plot)=="immbernie"] <- "Sanders (immig)"
rownames(df.plot)[rownames(df.plot)=="immtrump"] <- "Trump (immig)"
rownames(df.plot)[rownames(df.plot)=="cultrdem"] <- "DParty (cultural)"
rownames(df.plot)[rownames(df.plot)=="cultrrep"] <- "RParty (cultural)"
rownames(df.plot)[rownames(df.plot)=="cultrbernie"] <- "Sanders (cultural)"
rownames(df.plot)[rownames(df.plot)=="cultrtrump"] <- "Trump (cultural)"
rownames(df.plot)[rownames(df.plot)=="envdem"] <- "DParty (environ)"
rownames(df.plot)[rownames(df.plot)=="envrep"] <- "RParty (environ)"
rownames(df.plot)[rownames(df.plot)=="envbernie"] <- "Sanders (environ)"
rownames(df.plot)[rownames(df.plot)=="envtrump"] <- "Trump (environ)"
rownames(df.plot)[rownames(df.plot)=="discdem"] <- "DParty (discrim)"
rownames(df.plot)[rownames(df.plot)=="discrep"] <- "RParty (discrim)"
rownames(df.plot)[rownames(df.plot)=="discbernie"] <- "Sanders (discrim)"
rownames(df.plot)[rownames(df.plot)=="disctrump"] <- "Trump (discrim)"
rownames(df.plot)[rownames(df.plot)=="healthdem"] <- "DParty (health)"
rownames(df.plot)[rownames(df.plot)=="healthrep"] <- "RParty (health)"
rownames(df.plot)[rownames(df.plot)=="healthbernie"] <- "Sanders (health)"
rownames(df.plot)[rownames(df.plot)=="healthtrump"] <- "Trump (health)"
rownames(df.plot)[rownames(df.plot)=="intgdem"] <- "DParty (ininteg)"
rownames(df.plot)[rownames(df.plot)=="intgrep"] <- "RParty (ininteg)"
rownames(df.plot)[rownames(df.plot)=="intgbernie"] <- "Sanders (ininteg)"
rownames(df.plot)[rownames(df.plot)=="intgtrump"] <- "Trump (ininteg)"
#
#
# Table 1: Stimuli producing significant differences in mean ratings of Americanist and Transatlantic experts.
#
print(df.plot[df.plot$pvalue < 0.1,])
#
#
# Figure 1: Difference of means tests between US and European expert placements of stimuli.
#
pdf("fig1.pdf", height=7.5, width=6)
#
ggplot(df.plot, aes(order, pvalue, xend=order, yend=0)) +
   geom_segment(colour="grey50") +
   geom_point(size=2, aes(colour="blue2")) +
#  geom_bar(stat = "identity", show.legend = FALSE) +
  # Free the scales here
#  facet_wrap(~ party, scales = "free") +
  xlab("") +
  ylab("\np-value") +
  ggtitle("2017 CHES-US Pilot Study\n") +
  ylim(0,1) +
scale_color_viridis(discrete=TRUE, option="plasma", begin = 0.0, end = 0.9) +
  theme_bw() +
  theme(legend.position="none", plot.title = element_text(hjust = 0.45), axis.title.x = element_text(hjust = 0.45)) +
   theme(panel.grid.major.x = element_blank()) +
  scale_x_continuous(
    breaks = df.plot$order,
   labels = rownames(df.plot),
    expand = c(0.025,0)) +
  coord_flip() +
  geom_hline(aes(yintercept=0.1))
#
dev.off()
#
#
# Simulation tests for p-values
#
nstimuli <- nrow(US.expert.placements)
US.fitted.normals <- matrix(NA, nrow=nstimuli, ncol=2)
EU.fitted.normals <- matrix(NA, nrow=nstimuli, ncol=2)
#
for (i in 1:nstimuli){
	US.fitted.normals[i,] <- MASS::fitdistr(na.omit(US.expert.placements[i,]), "normal")$estimate
	EU.fitted.normals[i,] <- MASS::fitdistr(na.omit(EU.expert.placements[i,]), "normal")$estimate
}
#
US.fitted.normals.0.5 <- US.fitted.normals
US.fitted.normals.1.0 <- US.fitted.normals
US.fitted.normals.1.5 <- US.fitted.normals
US.fitted.normals.2.0 <- US.fitted.normals
#
for (i in 1:nstimuli){
	set.seed(1985+i)
	US.fitted.normals.0.5[,1] <- EU.fitted.normals[,1] + sample(c(-0.5,0.5), 1)
	US.fitted.normals.1.0[,1] <- EU.fitted.normals[,1] + sample(c(-1.0,1.0), 1)
	US.fitted.normals.1.5[,1] <- EU.fitted.normals[,1] + sample(c(-1.5,1.5), 1)
	US.fitted.normals.2.0[,1] <- EU.fitted.normals[,1] + sample(c(-2.0,2.0), 1)
	}
#
#
simulated.pvalues.0.5 <- matrix(NA, nrow=nstimuli, ncol=1000)
simulated.pvalues.1.0 <- matrix(NA, nrow=nstimuli, ncol=1000)
simulated.pvalues.1.5 <- matrix(NA, nrow=nstimuli, ncol=1000)
simulated.pvalues.2.0 <- matrix(NA, nrow=nstimuli, ncol=1000)
#
for (b in 1:1000){
	for (i in 1:nstimuli){
		set.seed(1985*b+i)
		US.temp <- round(rnorm(13, US.fitted.normals.0.5[i,1], US.fitted.normals.0.5[i,2]))
		EU.temp <- round(rnorm(25, EU.fitted.normals[i,1], EU.fitted.normals[i,2]))
		US.temp[US.temp < 1] <- 1
  		US.temp[US.temp > 10] <- 10
		EU.temp[EU.temp < 1] <- 1
  		EU.temp[EU.temp > 10] <- 10
		simulated.pvalues.0.5[i,b] <- t.test(US.temp, EU.temp)$p.value
	}}
#
for (b in 1:1000){
	for (i in 1:nstimuli){
		set.seed(1985*b+i)
		US.temp <- round(rnorm(13, US.fitted.normals.1.0[i,1], US.fitted.normals.1.0[i,2]))
		EU.temp <- round(rnorm(25, EU.fitted.normals[i,1], EU.fitted.normals[i,2]))
		US.temp[US.temp < 1] <- 1
  		US.temp[US.temp > 10] <- 10
		EU.temp[EU.temp < 1] <- 1
  		EU.temp[EU.temp > 10] <- 10
		simulated.pvalues.1.0[i,b] <- t.test(US.temp, EU.temp)$p.value
	}}
#
for (b in 1:1000){
	for (i in 1:nstimuli){
		set.seed(1985*b+i)
		US.temp <- round(rnorm(13, US.fitted.normals.1.5[i,1], US.fitted.normals.1.5[i,2]))
		EU.temp <- round(rnorm(25, EU.fitted.normals[i,1], EU.fitted.normals[i,2]))
		US.temp[US.temp < 1] <- 1
  		US.temp[US.temp > 10] <- 10
		EU.temp[EU.temp < 1] <- 1
  		EU.temp[EU.temp > 10] <- 10
		simulated.pvalues.1.5[i,b] <- t.test(US.temp, EU.temp)$p.value
	}}
#
for (b in 1:1000){
	for (i in 1:nstimuli){
		set.seed(1985*b+i)
		US.temp <- round(rnorm(13, US.fitted.normals.2.0[i,1], US.fitted.normals.2.0[i,2]))
		EU.temp <- round(rnorm(25, EU.fitted.normals[i,1], EU.fitted.normals[i,2]))
		US.temp[US.temp < 1] <- 1
  		US.temp[US.temp > 10] <- 10
		EU.temp[EU.temp < 1] <- 1
  		EU.temp[EU.temp > 10] <- 10
		simulated.pvalues.2.0[i,b] <- t.test(US.temp, EU.temp)$p.value
	}}
#
#
proptrue.0.5 <- apply(simulated.pvalues.0.5, 1, function(x){table(x < 0.1)["TRUE"] / length(x)})
proptrue.1.0 <- apply(simulated.pvalues.1.0, 1, function(x){table(x < 0.1)["TRUE"] / length(x)})
proptrue.1.5 <- apply(simulated.pvalues.1.5, 1, function(x){table(x < 0.1)["TRUE"] / length(x)})
proptrue.2.0 <- apply(simulated.pvalues.2.0, 1, function(x){table(x < 0.1)["TRUE"] / length(x)})
#
sim.pvalues.table <- cbind(proptrue.0.5, proptrue.1.0, proptrue.1.5, proptrue.2.0)
rownames(sim.pvalues.table) <- rownames(US.expert.placements)
#
#
#
#
# BASIC DEMOGRAPHICS
#
country <- as.factor(CHESUS$country)
#
leftrightself <- CHESUS$leftrightself
#
#
# VIGNETTE PLACEMENTS
#
partyAecon <- CHESUS$partyaecon
partyBecon <- CHESUS$partybecon
partyCecon <- CHESUS$partycecon
#
partyAsocial <- CHESUS$partyasocial
partyBsocial <- CHESUS$partybsocial
partyCsocial <- CHESUS$partycsocial
#
#
# US STIMULI PLACEMENTS
#
econ.US.dem <- CHESUS$econdem
econ.US.rep <- CHESUS$econrep
econ.US.sanders <- CHESUS$econbernie
econ.US.trump <- CHESUS$econtrump

social.US.dem <- CHESUS$socialdem
social.US.rep <- CHESUS$socialrep
social.US.sanders <- CHESUS$socialbernie
social.US.trump <- CHESUS$socialtrump
#
#
# EUROPEAN STIMULI PLACEMENTS
#
econ.FR.socialists <- CHESUS$econfrsoc
econ.FR.republicans <- CHESUS$econfrrep
econ.FR.left <- CHESUS$econfrleft
econ.FR.natfront <- CHESUS$econfrnatl
econ.FR.greens <- CHESUS$econfreco
#
econ.UK.labour <- CHESUS$econlab
econ.UK.conservatives <- CHESUS$econcons
econ.UK.libdems <- CHESUS$econlib
econ.UK.ukip <- CHESUS$econukip
econ.UK.greens <- CHESUS$econgreen
#
econ.DE.cdu <- CHESUS$econcdu
econ.DE.csu <- CHESUS$econcsu
econ.DE.spd <- CHESUS$econspd
econ.DE.fdp <- CHESUS$econfdp
econ.DE.grun <- CHESUS$econgrun
econ.DE.linke <- CHESUS$econlinke
econ.DE.afg <- CHESUS$econafg
#
social.FR.socialists <- CHESUS$socialsoc
social.FR.republicans <- CHESUS$socialrepub
social.FR.left <- CHESUS$socialleft
social.FR.natfront <- CHESUS$socialnatl
social.FR.greens <- CHESUS$socialeco
#
social.UK.labour <- CHESUS$sociallab
social.UK.conservatives <- CHESUS$socialcons
social.UK.libdems <- CHESUS$sociallib
social.UK.ukip <- CHESUS$socialukip
social.UK.greens <- CHESUS$socialgreen
#
social.DE.cdu <- CHESUS$socialcdu
social.DE.csu <- CHESUS$socialcsu
social.DE.spd <- CHESUS$socialspd
social.DE.fdp <- CHESUS$socialfdp
social.DE.grun <- CHESUS$socialgrun
social.DE.linke <- CHESUS$sociallinke
social.DE.afg <- CHESUS$socialafg
#
#
#
#
#  BAYESIAN ALDRICH-MCKELVEY SCALING
#    ******ECONOMIC ISSUES******
#
correct.ordering <- as.numeric((partyBecon < partyAecon) & (partyAecon < partyCecon))
correct.ordering[is.na(correct.ordering)] <- 1
#
placements <- data.frame(partyAecon, partyBecon, partyCecon,
	econ.US.dem, econ.US.rep, econ.US.sanders, econ.US.trump,
	econ.FR.socialists, econ.FR.republicans, econ.FR.left, econ.FR.natfront, econ.FR.greens,
	econ.UK.labour, econ.UK.conservatives, econ.UK.libdems, econ.UK.ukip, econ.UK.greens,
	econ.DE.cdu, econ.DE.csu, econ.DE.spd, econ.DE.fdp, econ.DE.grun, econ.DE.linke, econ.DE.afg)
placements[placements < 0 | placements > 10] <- NA
#
us.stimuli.names <- c("econ.US.dem", "econ.US.rep", "econ.US.sanders", "econ.US.trump")
vignette.names <- c("partyAecon", "partyBecon", "partyCecon")
#
# (1)
just.US.experts <- placements[country=="US", 1:7]
#
# (2)
vignettesonly <- placements
vignettesonly[country!="US",names(vignettesonly) %in% us.stimuli.names] <- NA
#
# (3)
novignettes <- placements
novignettes <- novignettes[,!names(novignettes) %in% vignette.names]
#
# (4)
fulldata <- placements
#
# (5)
justDR <- novignettes
justDR[country!="US",names(justDR) %in% c("econ.US.sanders", "econ.US.trump")] <- NA
justDR.economic <- justDR
#
# 5 kinds of datasets:
# (1) Just US experts
# (2) removing European placements of US stimuli (vignettes only)
# (3) removing vignettes (no vignettes)
# (4) keeping European placements of US stimuli (full data)
# (5) Just Dems and Reps as bridges
#
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# 1.) Just US experts
#
placementdat <- tempjunk <- data.frame(just.US.experts)
#
cutoff <- 4   # This section removes respondents who placed less than 4 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(tempjunk)) >= cutoff & apply(tempjunk, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- placementdat[index,]
#
TT <- T - 5
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
for (j in 1:3){
zhat[j] ~ dnorm(0,1)
}
zhat[4] ~ dnorm(0,1)T(-1.1,-1.0)
zhat[5] ~ dnorm(0,1)T(0.9,1.1)
for (j in 6:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(0,-1,1,-1,1,rnorm(q-5))
)}
#
BAM.sim <- jags.model(textConnection(modelstring),
    data = list('z' = z, 'q' = q, 'N' = N),
    inits = inits, n.chains = 2, n.adapt = 50000)
#
update(BAM.sim, n.iter=50000)
#
params <- c("zhat", "a", "b")
samps <- coda.samples(BAM.sim, params, n.iter=100000, thin=20)
#
zhatsamples.fullplacements <- do.call("rbind", samps[,grep(pattern="zhat[", colnames(samps[[1]]), fixed=TRUE)])
#
BAM1 <- colMeans(zhatsamples.fullplacements)
names(BAM1) <- colnames(T)
#
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# 2.) Vignettes only
#
placementdat <- tempjunk <- data.frame(vignettesonly)
#
cutoff <- 4   # This section removes respondents who placed less than 4 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(tempjunk)) >= cutoff & apply(tempjunk, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- placementdat[index,]
#
TT <- T - 5
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
for (j in 1:3){
zhat[j] ~ dnorm(0,1)
}
zhat[4] ~ dnorm(0,1)T(-1.1,-1.0)
zhat[5] ~ dnorm(0,1)T(0.9,1.1)
for (j in 6:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(0,-1,1,-1,1,rnorm(q-5))
)}
#
BAM.sim <- jags.model(textConnection(modelstring),
    data = list('z' = z, 'q' = q, 'N' = N),
    inits = inits, n.chains = 2, n.adapt = 50000)
#
update(BAM.sim, n.iter=50000)
#
params <- c("zhat", "a", "b")
samps <- coda.samples(BAM.sim, params, n.iter=100000, thin=20)
#
alphasamples.fullplacements <- do.call("rbind", samps[,grep(pattern="a[", colnames(samps[[1]]), fixed=TRUE)])
betasamples.fullplacements <- do.call("rbind", samps[,grep(pattern="b[", colnames(samps[[1]]), fixed=TRUE)])
#
beta.economic.FR <- data.frame(country="France", beta=as.vector(betasamples.fullplacements[,country[index]=="FR"]))
beta.economic.GE <- data.frame(country="Germany", beta=as.vector(betasamples.fullplacements[,country[index]=="GE"]))
beta.economic.UK <- data.frame(country="UK", beta=as.vector(betasamples.fullplacements[,country[index]=="UK"]))
beta.economic.US <- data.frame(country="US", beta=as.vector(betasamples.fullplacements[,country[index]=="US"]))
#
betas.economic <- rbind(beta.economic.FR, beta.economic.GE, beta.economic.UK, beta.economic.US)
#
aggregate(betas.economic$beta, list(betas.economic$country), mean)
t.test(betas.economic$beta[betas.economic$country=="US"], betas.economic$beta[betas.economic$country!="US"])
#
# Figure A2. Samples from posterior distributions of expert i values by country of expertise for
# the general economic left-right scale. Results are from Bayesian Aldrich-McKelvey scaling of the
# ‘vignettes only’ and ‘full data’ specifications, using the 36 experts who provided at least four valid
# stimuli placements. County-specific means shown.
#
ann.text <- data.frame(country=levels(betas.economic$country),
	  beta=0, y=0.5,
	  lab=c(mean(betas.economic$beta[betas.economic$country=="France"]),
			mean(betas.economic$beta[betas.economic$country=="Germany"]),
			mean(betas.economic$beta[betas.economic$country=="UK"]),
			mean(betas.economic$beta[betas.economic$country=="US"])))
#
pdf("figA2a.pdf", height=6.5, width=7)
#
ggplot(betas.economic, aes(x=beta, group=country, color=country, fill=country)) +
	geom_histogram(aes(y=..density..), position="identity", bins=50, alpha=0.4) +
	facet_wrap(~country) +
	theme(legend.position="none") +
	ggtitle("Stretch Parameter (Beta) by Country\nEconomic Issues - Vignettes Only") +
	ylab("") +
	xlab("Beta (sampled posterior values)") +
  geom_text(data = ann.text, aes(x=beta-2,y=y-0.005, group=country, size=0.75), label = paste("mu"), parse = TRUE) +
  geom_text(data = ann.text, aes(y=y, group=country, size=0.75), label = paste("=", round(ann.text$lab,2)), parse = FALSE)
#
dev.off()
#
zhatsamples.fullplacements <- do.call("rbind", samps[,grep(pattern="zhat[", colnames(samps[[1]]), fixed=TRUE)])
#
BAM2 <- colMeans(zhatsamples.fullplacements)
names(BAM2) <- colnames(T)
#
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# 3.) No vignettes
#
placementdat <- tempjunk <- data.frame(novignettes)
#
cutoff <- 4   # This section removes respondents who placed less than 4 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(tempjunk)) >= cutoff & apply(tempjunk, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- placementdat[index,]
#
TT <- T - 5
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
zhat[1] ~ dnorm(0,1)T(-1.1,-1.0)
zhat[2] ~ dnorm(0,1)T(0.9,1.1)
for (j in 3:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(-1,1,rnorm(q-2))
)}
#
BAM.sim <- jags.model(textConnection(modelstring),
    data = list('z' = z, 'q' = q, 'N' = N),
    inits = inits, n.chains = 2, n.adapt = 50000)
#
update(BAM.sim, n.iter=50000)
#
params <- c("zhat", "a", "b")
samps <- coda.samples(BAM.sim, params, n.iter=100000, thin=20)
#
#
#
zhatsamples.fullplacements <- do.call("rbind", samps[,grep(pattern="zhat[", colnames(samps[[1]]), fixed=TRUE)])
#
BAM3 <- colMeans(zhatsamples.fullplacements)
names(BAM3) <- colnames(T)
BAM3.economic <- BAM3
#
####################################################
####################################################
#
# 4.) Full data
#
placementdat <- tempjunk <- data.frame(fulldata)
#
cutoff <- 4   # This section removes respondents who placed less than 4 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(tempjunk)) >= cutoff & apply(tempjunk, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- placementdat[index,]
#
TT <- T - 5
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
for (j in 1:3){
zhat[j] ~ dnorm(0,1)
}
zhat[4] ~ dnorm(0,1)T(-1.1,-1.0)
zhat[5] ~ dnorm(0,1)T(0.9,1.1)
for (j in 6:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(0,-1,1,-1,1,rnorm(q-5))
)}
#
BAM.sim <- jags.model(textConnection(modelstring),
    data = list('z' = z, 'q' = q, 'N' = N),
    inits = inits, n.chains = 2, n.adapt = 50000)
#
update(BAM.sim, n.iter=50000)
#
params <- c("zhat", "a", "b")
samps <- coda.samples(BAM.sim, params, n.iter=100000, thin=20)
#
#
alphasamples.fullplacements <- do.call("rbind", samps[,grep(pattern="a[", colnames(samps[[1]]), fixed=TRUE)])
betasamples.fullplacements <- do.call("rbind", samps[,grep(pattern="b[", colnames(samps[[1]]), fixed=TRUE)])
#
beta.economic.FR <- data.frame(country="France", beta=as.vector(betasamples.fullplacements[,country[index]=="FR"]))
beta.economic.GE <- data.frame(country="Germany", beta=as.vector(betasamples.fullplacements[,country[index]=="GE"]))
beta.economic.UK <- data.frame(country="UK", beta=as.vector(betasamples.fullplacements[,country[index]=="UK"]))
beta.economic.US <- data.frame(country="US", beta=as.vector(betasamples.fullplacements[,country[index]=="US"]))
#
betas.economic <- rbind(beta.economic.FR, beta.economic.GE, beta.economic.UK, beta.economic.US)
#
t.test(betas.economic$beta[betas.economic$country=="US"], betas.economic$beta[betas.economic$country!="US"])
aggregate(betas.economic$beta, list(betas.economic$country), mean)
#
#
# Figure A2. Samples from posterior distributions of experti values by country of expertise for
# the general economic left-right scale. Results are from Bayesian Aldrich-McKelvey scaling of the
# ‘vignettes only’ and ‘full data’ specifications, using the 36 experts who provided at least four valid
# stimuli placements. County-specific means shown.
#
ann.text <- data.frame(country=levels(betas.economic$country),
	  beta=0, y=0.5,
	  lab=c(mean(betas.economic$beta[betas.economic$country=="France"]),
			mean(betas.economic$beta[betas.economic$country=="Germany"]),
			mean(betas.economic$beta[betas.economic$country=="UK"]),
			mean(betas.economic$beta[betas.economic$country=="US"])))
#
pdf("figA2b.pdf", height=6.5, width=7)
#
ggplot(betas.economic, aes(x=beta, group=country, color=country, fill=country)) +
	geom_histogram(aes(y=..density..), position="identity", bins=50, alpha=0.4) +
	facet_wrap(~country) +
	theme(legend.position="none") +
	ggtitle("Stretch Parameter (Beta) by Country\nEconomic Issues - Full Data") +
	ylab("") +
	xlab("Beta (sampled posterior values)") +
  geom_text(data = ann.text, aes(x=beta-2,y=y-0.005, group=country, size=0.75), label = paste("mu"), parse = TRUE) +
  geom_text(data = ann.text, aes(y=y, group=country, size=0.75), label = paste("=", round(ann.text$lab,2)), parse = FALSE)
#
#
dev.off()
#
zhatsamples.fullplacements <- do.call("rbind", samps[,grep(pattern="zhat[", colnames(samps[[1]]), fixed=TRUE)])
#
BAM4 <- colMeans(zhatsamples.fullplacements)
names(BAM4) <- colnames(T)
#
#
#
#
justUS <- c(BAM1, rep(NA, 17))
vignettesonly <- BAM2
novignettes <- c(rep(NA, 3), BAM3)
full <- BAM4
#
dat1 <- data.frame(Americanists=justUS, vignettes.only=vignettesonly, no.vignettes=novignettes, full.data=full, id=names(BAM2))
dat2 <- melt(dat1)
#
levels(dat2$variable) <- c("Americanists", "Vignettes only", "No vignettes", "Full data")
labs <- names(BAM2)
#
labs[labs=="partyAecon"] <- "Party A"
labs[labs=="partyBecon"] <- "Party B"
labs[labs=="partyCecon"] <- "Party C"
labs[labs=="econ.US.dem"] <- "D-Party (us)"
labs[labs=="econ.US.rep"] <- "R-Party (us)"
labs[labs=="econ.US.sanders"] <- "Sanders (us)"
labs[labs=="econ.US.trump"] <- "Trump (us)"
labs[labs=="econ.FR.socialists"] <- "Socialists (fr)"
labs[labs=="econ.FR.republicans"] <- "Repub (fr)"
labs[labs=="econ.FR.left"] <- "Left (fr)"
labs[labs=="econ.FR.natfront"] <- "Nat Front (fr)"
labs[labs=="econ.FR.greens"] <- "Greens (fr)"
labs[labs=="econ.UK.labour"] <- "Labour (uk)"
labs[labs=="econ.UK.conservatives"] <- "Conserv (uk)"
labs[labs=="econ.UK.libdems"] <- "LibDems (uk)"
labs[labs=="econ.UK.ukip"] <- "UKIP (uk)"
labs[labs=="econ.UK.greens"] <- "Greens (uk)"
labs[labs=="econ.DE.cdu"] <- "CDU (de)"
labs[labs=="econ.DE.csu"] <- "CSU (de)"
labs[labs=="econ.DE.spd"] <- "SPD (de)"
labs[labs=="econ.DE.fdp"] <- "FDP (de)"
labs[labs=="econ.DE.grun"] <- "Greens (de)"
labs[labs=="econ.DE.linke"] <- "Linke (de)"
labs[labs=="econ.DE.afg"] <- "AfD (de)"
#
labelsL <- paste(labs, "   ")
labelsR <- paste("   ", labs)
#
labelsL[8:24] <- ""
labelsR[1:7] <- ""
#
colors <- rep(viridis(24, begin=0.1, end=0.9, direction=1), 4)
dat2$colors <- rep(viridis(24, begin=0.1, end=0.9, direction=1), 4)
#
# Figure 2. Bayesian Aldrich-McKelvey ideal point estimates of general economic positions under alternative data specifications
#
pdf("fig2.pdf", height=9.8, width=9.75)
#
ggplot(dat2) +
  geom_line(aes(x = as.factor(variable), y = value, group = id, color=colors), size = 2) +
  geom_point(aes(x = as.factor(variable), y = value, color=colors), size = 5) +
  theme_minimal(base_size = 18) +
  scale_color_manual(values=c(dat2$colors[1:24])) +
  #scale_color_manual(values=rep(viridis(24, begin=0.1, end=0.9, direction=1), 4)) +
  #scale_color_viridis(discrete=TRUE, option="C", end=0.9, direction=1)   +
  #scale_color_distiller(palette = "Spectral") +
xlab("") +
geom_text(data = subset(dat2, variable=="Americanists"),
		aes(x = as.factor(variable), y = value, color = colors, label = labelsL),
        size = 5, hjust = 1, vjust=c(0,0,0,0,0,0,0,0,
									 0,0,0,0,0,0,0,0,
									 0,0,0,0,0,0,0,0)) +
geom_text(data = subset(dat2, variable=="Full data"),
        aes(x = as.factor(variable), y = value, color = colors, label = labelsR),
        size = 5, hjust = 0, vjust=c(0,0,0,0,0,0,0,0.5,
									 1,0.5,-0.5,0,0,0,0,0,
									 0.5,0,0,0,0,1,-0.5,0)) +
  theme(legend.position = "none",
#      panel.grid.major.y = element_blank(),
      panel.grid.minor.y = element_blank(),
#      panel.grid.major.x = element_blank(),
      axis.ticks.y = element_blank(),
      axis.ticks.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.y = element_blank(),
      plot.title = element_text(size = 18, hjust=0.5)) +
      ggtitle("Economic left-right estimates")
#
dev.off()
#
#
#
#
#  BAYESIAN ALDRICH-MCKELVEY SCALING
#    ******SOCIAL ISSUES******
#
correct.ordering <- as.numeric((partyBsocial < partyAsocial) & (partyAsocial < partyCsocial))
correct.ordering[is.na(correct.ordering)] <- 1
#
placements <- data.frame(partyAsocial, partyBsocial, partyCsocial,
	social.US.dem, social.US.rep, social.US.sanders, social.US.trump,
	social.FR.socialists, social.FR.republicans, social.FR.left, social.FR.natfront, social.FR.greens,
	social.UK.labour, social.UK.conservatives, social.UK.libdems, social.UK.ukip, social.UK.greens,
	social.DE.cdu, social.DE.csu, social.DE.spd, social.DE.fdp, social.DE.grun, social.DE.linke, social.DE.afg)
placements[placements < 0 | placements > 10] <- NA
#
us.stimuli.names <- c("social.US.dem", "social.US.rep", "social.US.sanders", "social.US.trump")
vignette.names <- c("partyAsocial", "partyBsocial", "partyCsocial")
#
# (1)
just.US.experts <- placements[country=="US", 1:7]
#
# (2)
vignettesonly <- placements
vignettesonly[country!="US",names(vignettesonly) %in% us.stimuli.names] <- NA
#
# (3)
novignettes <- placements
novignettes <- novignettes[,!names(novignettes) %in% vignette.names]
#
# (4)
fulldata <- placements
#
# (5)
justDR <- novignettes
justDR[country!="US",names(justDR) %in% c("social.US.sanders", "social.US.trump")] <- NA
justDR.social <- justDR
#
# 5 kinds of datasets:
# (1) Just US experts
# (2) removing European placements of US stimuli (vignettes only)
# (3) removing vignettes (no vignettes)
# (4) keeping European placements of US stimuli (full data)
# (5) Just Dems and Reps as bridges
#
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# 1.) Just US experts
#
placementdat <- tempjunk <- data.frame(just.US.experts)
#
cutoff <- 4   # This section removes respondents who placed less than 4 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(tempjunk)) >= cutoff & apply(tempjunk, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- placementdat[index,]
#
TT <- T - 5
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
for (j in 1:3){
zhat[j] ~ dnorm(0,1)
}
zhat[4] ~ dnorm(0,1)T(-1.1,-1.0)
zhat[5] ~ dnorm(0,1)T(0.9,1.1)
for (j in 6:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(0,-1,1,-1,1,rnorm(q-5))
)}
#
BAM.sim <- jags.model(textConnection(modelstring),
    data = list('z' = z, 'q' = q, 'N' = N),
    inits = inits, n.chains = 2, n.adapt = 50000)
#
update(BAM.sim, n.iter=50000)
#
params <- c("zhat", "a", "b")
samps <- coda.samples(BAM.sim, params, n.iter=100000, thin=20)
#
zhatsamples.fullplacements <- do.call("rbind", samps[,grep(pattern="zhat[", colnames(samps[[1]]), fixed=TRUE)])
#
BAM1 <- colMeans(zhatsamples.fullplacements)
names(BAM1) <- colnames(T)
#
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# 2.) Vignettes only
#
placementdat <- tempjunk <- data.frame(vignettesonly)
#
cutoff <- 4   # This section removes respondents who placed less than 4 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(tempjunk)) >= cutoff & apply(tempjunk, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- placementdat[index,]
#
TT <- T - 5
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
for (j in 1:3){
zhat[j] ~ dnorm(0,1)
}
zhat[4] ~ dnorm(0,1)T(-1.1,-1.0)
zhat[5] ~ dnorm(0,1)T(0.9,1.1)
for (j in 6:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(0,-1,1,-1,1,rnorm(q-5))
)}
#
BAM.sim <- jags.model(textConnection(modelstring),
    data = list('z' = z, 'q' = q, 'N' = N),
    inits = inits, n.chains = 2, n.adapt = 50000)
#
update(BAM.sim, n.iter=50000)
#
params <- c("zhat", "a", "b")
samps <- coda.samples(BAM.sim, params, n.iter=100000, thin=20)
#
alphasamples.fullplacements <- do.call("rbind", samps[,grep(pattern="a[", colnames(samps[[1]]), fixed=TRUE)])
betasamples.fullplacements <- do.call("rbind", samps[,grep(pattern="b[", colnames(samps[[1]]), fixed=TRUE)])
#
beta.social.FR <- data.frame(country="France", beta=as.vector(betasamples.fullplacements[,country[index]=="FR"]))
beta.social.GE <- data.frame(country="Germany", beta=as.vector(betasamples.fullplacements[,country[index]=="GE"]))
beta.social.UK <- data.frame(country="UK", beta=as.vector(betasamples.fullplacements[,country[index]=="UK"]))
beta.social.US <- data.frame(country="US", beta=as.vector(betasamples.fullplacements[,country[index]=="US"]))
#
betas.social <- rbind(beta.social.FR, beta.social.GE, beta.social.UK, beta.social.US)
#
aggregate(betas.social$beta, list(betas.social$country), mean)
t.test(betas.social$beta[betas.social$country=="US"], betas.social$beta[betas.social$country!="US"])
#
# Figure A3. Samples from posterior distributions of expert i values by country of expertise for the
# general social left-right scale. Results are from Bayesian Aldrich-McKelvey scaling of the ‘vignettes
# only’ and ‘full data’ specifications, using the 35 experts who provided at least four valid stimuli
# placements. County-specific means shown.
#
ann.text <- data.frame(country=levels(betas.social$country),
	  beta=0, y=0.5,
	  lab=c(mean(betas.social$beta[betas.social$country=="France"]),
			mean(betas.social$beta[betas.social$country=="Germany"]),
			mean(betas.social$beta[betas.social$country=="UK"]),
			mean(betas.social$beta[betas.social$country=="US"])))
#
pdf("figA3a.pdf", height=6.5, width=7)
#
ggplot(betas.social, aes(x=beta, group=country, color=country, fill=country)) +
	geom_histogram(aes(y=..density..), position="identity", bins=50, alpha=0.4) +
	facet_wrap(~country) +
	theme(legend.position="none") +
	ggtitle("Stretch Parameter (Beta) by Country\nSocial Issues - Vignettes Only") +
	ylab("") +
	xlab("Beta (sampled posterior values)") +
  geom_text(data = ann.text, aes(x=beta-2,y=y-0.005, group=country, size=0.75), label = paste("mu"), parse = TRUE) +
  geom_text(data = ann.text, aes(y=y, group=country, size=0.75), label = paste("=", round(ann.text$lab,2)), parse = FALSE)
#
dev.off()
#
#
zhatsamples.fullplacements <- do.call("rbind", samps[,grep(pattern="zhat[", colnames(samps[[1]]), fixed=TRUE)])
#
BAM2 <- colMeans(zhatsamples.fullplacements)
names(BAM2) <- colnames(T)
#
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# 3.) No vignettes
#
placementdat <- tempjunk <- data.frame(novignettes)
#
cutoff <- 4   # This section removes respondents who placed less than 4 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(tempjunk)) >= cutoff & apply(tempjunk, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- placementdat[index,]
#
TT <- T - 5
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
zhat[1] ~ dnorm(0,1)T(-1.1,-1.0)
zhat[2] ~ dnorm(0,1)T(0.9,1.1)
for (j in 3:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(-1,1,rnorm(q-2))
)}
#
BAM.sim <- jags.model(textConnection(modelstring),
    data = list('z' = z, 'q' = q, 'N' = N),
    inits = inits, n.chains = 2, n.adapt = 50000)
#
update(BAM.sim, n.iter=50000)
#
params <- c("zhat", "a", "b")
samps <- coda.samples(BAM.sim, params, n.iter=100000, thin=20)
#
zhatsamples.fullplacements <- do.call("rbind", samps[,grep(pattern="zhat[", colnames(samps[[1]]), fixed=TRUE)])
#
BAM3 <- colMeans(zhatsamples.fullplacements)
names(BAM3) <- colnames(T)
BAM3.social <- BAM3
#
####################################################
####################################################
#
# 4.) Full data
#
placementdat <- tempjunk <- data.frame(fulldata)
#
cutoff <- 4   # This section removes respondents who placed less than 4 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(tempjunk)) >= cutoff & apply(tempjunk, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- placementdat[index,]
#
TT <- T - 5
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
for (j in 1:3){
zhat[j] ~ dnorm(0,1)
}
zhat[4] ~ dnorm(0,1)T(-1.1,-1.0)
zhat[5] ~ dnorm(0,1)T(0.9,1.1)
for (j in 6:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(0,-1,1,-1,1,rnorm(q-5))
)}
#
BAM.sim <- jags.model(textConnection(modelstring),
    data = list('z' = z, 'q' = q, 'N' = N),
    inits = inits, n.chains = 2, n.adapt = 50000)
#
update(BAM.sim, n.iter=50000)
#
params <- c("zhat", "a", "b")
samps <- coda.samples(BAM.sim, params, n.iter=100000, thin=20)
#
zhatsamples.fullplacements <- do.call("rbind", samps[,grep(pattern="zhat[", colnames(samps[[1]]), fixed=TRUE)])
#
BAM4 <- colMeans(zhatsamples.fullplacements)
names(BAM4) <- colnames(T)
#
alphasamples.fullplacements <- do.call("rbind", samps[,grep(pattern="a[", colnames(samps[[1]]), fixed=TRUE)])
betasamples.fullplacements <- do.call("rbind", samps[,grep(pattern="b[", colnames(samps[[1]]), fixed=TRUE)])
#
beta.social.FR <- data.frame(country="France", beta=as.vector(betasamples.fullplacements[,country[index]=="FR"]))
beta.social.GE <- data.frame(country="Germany", beta=as.vector(betasamples.fullplacements[,country[index]=="GE"]))
beta.social.UK <- data.frame(country="UK", beta=as.vector(betasamples.fullplacements[,country[index]=="UK"]))
beta.social.US <- data.frame(country="US", beta=as.vector(betasamples.fullplacements[,country[index]=="US"]))
#
betas.social <- rbind(beta.social.FR, beta.social.GE, beta.social.UK, beta.social.US)
#
aggregate(betas.social$beta, list(betas.social$country), mean)
t.test(betas.social$beta[betas.social$country=="US"], betas.social$beta[betas.social$country!="US"])
#
# Figure A3. Samples from posterior distributions of expert i values by country of expertise for the
# general social left-right scale. Results are from Bayesian Aldrich-McKelvey scaling of the ‘vignettes
# only’ and ‘full data’ specifications, using the 35 experts who provided at least four valid stimuli
# placements. County-specific means shown.
#
ann.text <- data.frame(country=levels(betas.social$country),
	  beta=0, y=0.5,
	  lab=c(mean(betas.social$beta[betas.social$country=="France"]),
			mean(betas.social$beta[betas.social$country=="Germany"]),
			mean(betas.social$beta[betas.social$country=="UK"]),
			mean(betas.social$beta[betas.social$country=="US"])))
#
pdf("figA3b.pdf", height=6.5, width=7)
#
ggplot(betas.social, aes(x=beta, group=country, color=country, fill=country)) +
	geom_histogram(aes(y=..density..), position="identity", bins=50, alpha=0.4) +
	facet_wrap(~country) +
	theme(legend.position="none") +
	ggtitle("Stretch Parameter (Beta) by Country\nSocial Issues - Full Data") +
	ylab("") +
	xlab("Beta (sampled posterior values)") +
  geom_text(data = ann.text, aes(x=beta-2,y=y-0.005, group=country, size=0.75), label = paste("mu"), parse = TRUE) +
  geom_text(data = ann.text, aes(y=y, group=country, size=0.75), label = paste("=", round(ann.text$lab,2)), parse = FALSE)
#
#
dev.off()
#
#
#
#
justUS <- c(BAM1, rep(NA, 17))
vignettesonly <- BAM2
novignettes <- c(rep(NA, 3), BAM3)
full <- BAM4
#
dat1 <- data.frame(Americanists=justUS, vignettes.only=vignettesonly, no.vignettes=novignettes, full.data=full, id=names(BAM2))
dat2 <- melt(dat1)
#
levels(dat2$variable) <- c("Americanists", "Vignettes only", "No vignettes", "Full data")
labs <- names(BAM2)
#
labs[labs=="partyAsocial"] <- "Party A"
labs[labs=="partyBsocial"] <- "Party B"
labs[labs=="partyCsocial"] <- "Party C"
labs[labs=="social.US.dem"] <- "D-Party (us)"
labs[labs=="social.US.rep"] <- "R-Party (us)"
labs[labs=="social.US.sanders"] <- "Sanders (us)"
labs[labs=="social.US.trump"] <- "Trump (us)"
labs[labs=="social.FR.socialists"] <- "Socialists (fr)"
labs[labs=="social.FR.republicans"] <- "Repub (fr)"
labs[labs=="social.FR.left"] <- "Left (fr)"
labs[labs=="social.FR.natfront"] <- "Nat Front (fr)"
labs[labs=="social.FR.greens"] <- "Greens (fr)"
labs[labs=="social.UK.labour"] <- "Labour (uk)"
labs[labs=="social.UK.conservatives"] <- "Conserv (uk)"
labs[labs=="social.UK.libdems"] <- "LibDems (uk)"
labs[labs=="social.UK.ukip"] <- "UKIP (uk)"
labs[labs=="social.UK.greens"] <- "Greens (uk)"
labs[labs=="social.DE.cdu"] <- "CDU (de)"
labs[labs=="social.DE.csu"] <- "CSU (de)"
labs[labs=="social.DE.spd"] <- "SPD (de)"
labs[labs=="social.DE.fdp"] <- "FDP (de)"
labs[labs=="social.DE.grun"] <- "Greens (de)"
labs[labs=="social.DE.linke"] <- "Linke (de)"
labs[labs=="social.DE.afg"] <- "AfD (de)"
#
labelsL <- paste(labs, "   ")
labelsR <- paste("   ", labs)
#
labelsL[8:24] <- ""
labelsR[1:7] <- ""
#
colors <- rep(viridis(24, begin=0.1, end=0.9, direction=1), 4)
dat2$colors <- rep(viridis(24, begin=0.1, end=0.9, direction=1), 4)
#
# Figure 3. Bayesian Aldrich-McKelvey ideal point estimates of general social/cultural positions
# under alternative data specifications.
#
pdf("fig3.pdf", height=9.8, width=9.75)
#
ggplot(dat2) +
  geom_line(aes(x = as.factor(variable), y = value, group = id, color=colors), size = 2) +
  geom_point(aes(x = as.factor(variable), y = value, color = colors), size = 5) +
  theme_minimal(base_size = 18) +
  scale_color_manual(values=c(dat2$colors[1:24])) +
  #scale_color_viridis(discrete=TRUE, option="C", end=0.9, direction=1)   +
  #scale_color_distiller(palette = "Spectral") +
xlab("") +
geom_text(data = subset(dat2, variable=="Americanists"),
		aes(x = as.factor(variable), y = value, color = colors, label = labelsL),
        size = 4.5, hjust = 1, vjust=c(-0.5,0,0.5,0.65,0,1.25,0,0,
									 0,0,0,0,0,0,0,0,
									 0,0,0,0,0,0,0,0)) +
geom_text(data = subset(dat2, variable=="Full data"),
        aes(x = as.factor(variable), y = value, color = colors, label = labelsR),
        size = 4.5, hjust = 0, vjust=c(0,0,0,0,0,0,0,0.75,
									 0,0,0,0.25,-0.5,0,0.5,0,
									 0.75,0,0,0,-0.5,-0.75,0.5,0)) +
  theme(legend.position = "none",
#      panel.grid.major.y = element_blank(),
      panel.grid.minor.y = element_blank(),
#      panel.grid.major.x = element_blank(),
      axis.ticks.y = element_blank(),
      axis.ticks.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.y = element_blank(),
      plot.title = element_text(size = 18, hjust=0.5)) +
      ggtitle("Social left-right estimates")
#
dev.off()
#
#
#
#
# INDIVIDUAL ISSUE SCALE PLACEMENTS
#
immigration <- CHESUS[,grep("imm", names(CHESUS))]
multiculturalism <- CHESUS[,grep("cultr", names(CHESUS))]
environment <- CHESUS[,grep("env", names(CHESUS))]
positivediscrimination <- CHESUS[,grep("disc", names(CHESUS))]
health <- CHESUS[,grep("health", names(CHESUS))]
internationalintegration <- CHESUS[,grep("intg", names(CHESUS))]
#
#
allissuescales <- cbind(immigration, multiculturalism,
	environment,positivediscrimination, health, internationalintegration)
#
#
# BAYESIAN ALDRICH-MCKELVEY SCALING FUNCTION
#
BAMfunction <- function(placementdat, cutoff=4){
	tempjunk <- placementdat
	index <- rowSums(!is.na(tempjunk)) >= cutoff & apply(tempjunk, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
	T <- placementdat[index,]
	TT <- T - 5
	N <- nrow(TT)
	q <- ncol(TT)
	z <- TT
	#
	modelstring <-
	"model{
		for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
		}}
		##priors on a and b
		for(i in 1:N){
			a[i] ~ dunif(-100,100)
			b[i] ~ dunif(-100,100)
		}
		##priors on variance
		for(j in 1:q){
		tauj[j] ~ dgamma(1,.1)
		}
		for(i in 1:N){
			taui[i] ~ dgamma(ga,gb)
		}
		ga ~ dgamma(1,.1)
		gb ~ dgamma(1,.1)
		##priors on zhat
		for (j in 1:q){
		zhat[j] ~ dnorm(0,1)
		}
		sigmai.sq <- 1/taui
		sigmaj.sq <- 1/tauj
		}"
#
	inits <- function(){
		list(
		a=runif(N,-2,2),
		b=runif(N,0.5,2),
		zhat=c(-1,1,((colMeans(z[,3:q], na.rm=TRUE) + rnorm(q-2)) / 10))
		)}
#
	BAM.sim <- jags.model(textConnection(modelstring),
    	data = list('z' = z, 'q' = q, 'N' = N),
    	inits = inits, n.chains = 2, n.adapt = 50000)
	update(BAM.sim, n.iter=50000)
	params <- c("zhat")
	samps <- coda.samples(BAM.sim, params, n.iter=100000, thin=20)
#
	samps.identified <- do.call("rbind", samps)
	samps.identified <- t(apply(samps.identified, 1, scale))
#
	zhatsamps <- samps.identified
#
	mean <- colMeans(zhatsamps)
	lb <- apply(zhatsamps, 2, function(x){quantile(x, probs=0.05)})
	ub <- apply(zhatsamps, 2, function(x){quantile(x, probs=0.95)})
	out <- data.frame(issue=mean,lb,ub)
	rownames(out) <- colnames(z)
#
return(out)
}
#
#
immigration.means <- BAMfunction(immigration)
multiculturalism.means <- BAMfunction(multiculturalism)
environment.means <- BAMfunction(environment)
positivediscrimination.means <- BAMfunction(positivediscrimination)
health.means <- BAMfunction(health)
internationalintegration.means <- BAMfunction(internationalintegration)
#
partynames <- c("D-Party (us)","R-Party (us)","Sanders (us)","Trump (us)",
"Socialists (fr)","Repub (fr)","Left (fr)","Nat Front (fr)","Greens (fr)",
"Labour (uk)","Conserv (uk)","LibDems (uk)","UKIP (uk)","Greens (uk)",
"CDU (de)","CSU (de)","SPD (de)","FDP (de)","Greens (de)","Linke (de)","AfD (de)")
#
A <- data.frame(partynames,
		Immigration=immigration.means[,1],
		Multiculturalism=multiculturalism.means[,1],
		Environment=environment.means[,1],
		Positive.Discrimination=positivediscrimination.means[,1],
		Health=health.means[,1],
		International.Integration=internationalintegration.means[,1])
#
maindf.tmp <- melt(A, id="partynames")
names(maindf.tmp) <- c("partynames", "issue", "mean")
#
A.lb <- data.frame(partynames,
		Immigration=immigration.means[,2],
		Multiculturalism=multiculturalism.means[,2],
		Environment=environment.means[,2],
		Positive.Discrimination=positivediscrimination.means[,2],
		Health=health.means[,2],
		International.Integration=internationalintegration.means[,2])
#
A.ub <- data.frame(partynames,
		Immigration=immigration.means[,3],
		Multiculturalism=multiculturalism.means[,3],
		Environment=environment.means[,3],
		Positive.Discrimination=positivediscrimination.means[,3],
		Health=health.means[,3],
		International.Integration=internationalintegration.means[,3])
#
maindf <- data.frame(maindf.tmp, melt(A.lb, id="partynames")[,3], melt(A.ub, id="partynames")[,3])
names(maindf) <- c("partynames", "issue", "idealpoint", "lb", "ub")
#
# subset by issue
intg = subset(maindf, issue=="International.Integration")
imm = subset(maindf, issue=="Immigration")
multclt = subset(maindf, issue=="Multiculturalism")
envnmt = subset(maindf, issue=="Environment")
pd = subset(maindf, issue=="Positive.Discrimination")
health = subset(maindf, issue=="Health")
#
#
intgPlot = qplot(x = idealpoint, y = reorder(partynames, idealpoint), data = intg, geom = "point") +
ylab("") +
xlab("Estimated ideal point") +
  geom_errorbarh(aes(xmin = lb, xmax = ub), height=0, size=0.3) +
  labs(title = "International Integration") +
  theme(plot.title = element_text(hjust = 0.5))
#
immPlot = qplot(x = idealpoint, y = reorder(partynames, idealpoint), data = imm, geom = "point") +
ylab("") +
xlab("Estimated ideal point") +
  geom_errorbarh(aes(xmin = lb, xmax = ub), height=0, size=0.3) +
  labs(title = "Immigration") +
  theme(plot.title = element_text(hjust = 0.5))
#
multcltPlot = qplot(x = idealpoint, y = reorder(partynames, idealpoint), data = multclt, geom = "point") +
ylab("") +
xlab("Estimated ideal point") +
  geom_errorbarh(aes(xmin = lb, xmax = ub), height=0, size=0.3) +
  labs(title = "Multiculturalism") +
  theme(plot.title = element_text(hjust = 0.5))
#
envnmtPlot = qplot(x = idealpoint, y = reorder(partynames, idealpoint), data = envnmt, geom = "point") +
ylab("") +
xlab("Estimated ideal point") +
  geom_errorbarh(aes(xmin = lb, xmax = ub), height=0, size=0.3) +
  labs(title = "Environment") +
  theme(plot.title = element_text(hjust = 0.5))
#
pdPlot = qplot(x = idealpoint, y = reorder(partynames, idealpoint), data = pd, geom = "point") +
ylab("") +
xlab("Estimated ideal point") +
  geom_errorbarh(aes(xmin = lb, xmax = ub), height=0, size=0.3) +
  labs(title = "Positive Discrimination") +
  theme(plot.title = element_text(hjust = 0.5))
#
healthPlot = qplot(x = idealpoint, y = reorder(partynames, idealpoint), data = health, geom = "point") +
  ylab("") +
  xlab("Estimated ideal point") +
  geom_errorbarh(aes(xmin = lb, xmax = ub), height=0, size=0.3) +
  labs(title = "Health") +
  theme(plot.title = element_text(hjust = 0.5))
#
# Figure 4. Policy-specific stimuli estimates.
#
pdf("fig4.pdf", height=10, width=7)
grid.arrange(intgPlot, immPlot, multcltPlot, envnmtPlot, pdPlot, healthPlot, ncol=2, nrow=3)
dev.off()
#
#
#
#
#
# justDR: Economic Issues
#
# 5.) justDR
#
placementdat <- tempjunk <- data.frame(justDR.economic)
#
cutoff <- 4   # This section removes respondents who placed less than 4 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(tempjunk)) >= cutoff & apply(tempjunk, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- placementdat[index,]
#
TT <- T - 5
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
zhat[1] ~ dnorm(0,1)T(-1.1,-1.0)
zhat[2] ~ dnorm(0,1)T(0.9,1.1)
for (j in 3:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(-1,1,rnorm(q-2))
)}
#
BAM.sim <- jags.model(textConnection(modelstring),
    data = list('z' = z, 'q' = q, 'N' = N),
    inits = inits, n.chains = 1, n.adapt = 50000)
#
update(BAM.sim, n.iter=50000)
#
params <- c("zhat", "a", "b")
samps <- coda.samples(BAM.sim, params, n.iter=100000, thin=20)
#
zhatsamples.fullplacements <- do.call("rbind", samps[,grep(pattern="zhat[", colnames(samps[[1]]), fixed=TRUE)])
#
BAM5 <- colMeans(zhatsamples.fullplacements)
names(BAM5) <- colnames(T)
#
#
#
novignettes <- c(rep(NA, 3), BAM3.economic)
onlyDR <- c(rep(NA, 3), BAM5)
#
dat1 <- data.frame(no.vignettes=novignettes, onlyDR=onlyDR, id=names(novignettes))
dat2 <- melt(dat1)
dat2 <- na.omit(dat2)
#
levels(dat2$variable) <- c("No vignettes", "Only DR")
labs <- names(novignettes)[-c(1,2,3)]
#
labs[labs=="econ.US.dem"] <- "D-Party (us)"
labs[labs=="econ.US.rep"] <- "R-Party (us)"
labs[labs=="econ.US.sanders"] <- "Sanders (us)"
labs[labs=="econ.US.trump"] <- "Trump (us)"
labs[labs=="econ.FR.socialists"] <- "Socialists (fr)"
labs[labs=="econ.FR.republicans"] <- "Repub (fr)"
labs[labs=="econ.FR.left"] <- "Left (fr)"
labs[labs=="econ.FR.natfront"] <- "Nat Front (fr)"
labs[labs=="econ.FR.greens"] <- "Greens (fr)"
labs[labs=="econ.UK.labour"] <- "Labour (uk)"
labs[labs=="econ.UK.conservatives"] <- "Conserv (uk)"
labs[labs=="econ.UK.libdems"] <- "LibDems (uk)"
labs[labs=="econ.UK.ukip"] <- "UKIP (uk)"
labs[labs=="econ.UK.greens"] <- "Greens (uk)"
labs[labs=="econ.DE.cdu"] <- "CDU (de)"
labs[labs=="econ.DE.csu"] <- "CSU (de)"
labs[labs=="econ.DE.spd"] <- "SPD (de)"
labs[labs=="econ.DE.fdp"] <- "FDP (de)"
labs[labs=="econ.DE.grun"] <- "Greens (de)"
labs[labs=="econ.DE.linke"] <- "Linke (de)"
labs[labs=="econ.DE.afg"] <- "AfD (de)"
#
labelsL <- paste(labs, "   ")
labelsR <- paste("   ", labs)
#
labelsL[12:21] <- ""
labelsR[1:11] <- ""
#
colors <- rep(viridis(21, begin=0.1, end=0.9, direction=1), 2)
dat2$colors <- rep(viridis(21, begin=0.1, end=0.9, direction=1), 2)
#
# Figure A1a. Results from alternate specifications of the joint Bayesian Aldrich-McKelvey model.
#
pdf("figA1a.pdf", height=9, width=8)
#
ggplot(dat2) +
  geom_line(aes(x = as.factor(variable), y = value, group = id, color=colors), size = 2) +
  geom_point(aes(x = as.factor(variable), y = value, color=colors), size = 5) +
  theme_minimal(base_size = 18) +
  scale_color_manual(values=c(dat2$colors[1:21])) +
  #scale_color_manual(values=rep(viridis(24, begin=0.1, end=0.9, direction=1), 4)) +
  #scale_color_viridis(discrete=TRUE, option="C", end=0.9, direction=1)   +
  #scale_color_distiller(palette = "Spectral") +
xlab("") +
geom_text(data = subset(dat2, variable=="No vignettes"),
		aes(x = as.factor(variable), y = value, color = colors, label = labelsL),
        size = 5, hjust = 1, vjust=c(0,0,0,0,0,
									 0,0,0,0,0,0,0,0,
									 0,0,0,0,0,0,0,0)) +
geom_text(data = subset(dat2, variable=="Only DR"),
        aes(x = as.factor(variable), y = value, color = colors, label = labelsR),
        size = 5, hjust = 0, vjust=c(0,0,0,0,0.5,
									 1,0.5,-0.5,0,0,0,0,0,
									 0.5,0,0,0,0,1,-0.5,0)) +
  theme(legend.position = "none",
#      panel.grid.major.y = element_blank(),
      panel.grid.minor.y = element_blank(),
#      panel.grid.major.x = element_blank(),
      axis.ticks.y = element_blank(),
      axis.ticks.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.y = element_blank(),
      plot.title = element_text(size = 18, hjust=0.5)) +
      ggtitle("Economic left-right estimates")
#
dev.off()
#
#
#
#
# justDR: Social Issues
#
# 5.) justDR
#
placementdat <- tempjunk <- data.frame(justDR.social)
#
cutoff <- 4   # This section removes respondents who placed less than 4 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(tempjunk)) >= cutoff & apply(tempjunk, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- placementdat[index,]
#
TT <- T - 5
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
zhat[1] ~ dnorm(0,1)T(-1.1,-1.0)
zhat[2] ~ dnorm(0,1)T(0.9,1.1)
for (j in 3:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(-1,1,rnorm(q-2))
)}
#
BAM.sim <- jags.model(textConnection(modelstring),
    data = list('z' = z, 'q' = q, 'N' = N),
    inits = inits, n.chains = 1, n.adapt = 50000)
#
update(BAM.sim, n.iter=50000)
#
params <- c("zhat", "a", "b")
samps <- coda.samples(BAM.sim, params, n.iter=100000, thin=20)
#
zhatsamples.fullplacements <- do.call("rbind", samps[,grep(pattern="zhat[", colnames(samps[[1]]), fixed=TRUE)])
#
BAM5 <- colMeans(zhatsamples.fullplacements)
names(BAM5) <- colnames(T)
#
#
#
novignettes <- c(rep(NA, 3), BAM3.social)
onlyDR <- c(rep(NA, 3), BAM5)
#
dat1 <- data.frame(no.vignettes=novignettes, onlyDR=onlyDR, id=names(novignettes))
dat2 <- melt(dat1)
dat2 <- na.omit(dat2)
#
levels(dat2$variable) <- c("No vignettes", "Only DR")
labs <- names(novignettes)[-c(1,2,3)]
#
labs[labs=="social.US.dem"] <- "D-Party (us)"
labs[labs=="social.US.rep"] <- "R-Party (us)"
labs[labs=="social.US.sanders"] <- "Sanders (us)"
labs[labs=="social.US.trump"] <- "Trump (us)"
labs[labs=="social.FR.socialists"] <- "Socialists (fr)"
labs[labs=="social.FR.republicans"] <- "Repub (fr)"
labs[labs=="social.FR.left"] <- "Left (fr)"
labs[labs=="social.FR.natfront"] <- "Nat Front (fr)"
labs[labs=="social.FR.greens"] <- "Greens (fr)"
labs[labs=="social.UK.labour"] <- "Labour (uk)"
labs[labs=="social.UK.conservatives"] <- "Conserv (uk)"
labs[labs=="social.UK.libdems"] <- "LibDems (uk)"
labs[labs=="social.UK.ukip"] <- "UKIP (uk)"
labs[labs=="social.UK.greens"] <- "Greens (uk)"
labs[labs=="social.DE.cdu"] <- "CDU (de)"
labs[labs=="social.DE.csu"] <- "CSU (de)"
labs[labs=="social.DE.spd"] <- "SPD (de)"
labs[labs=="social.DE.fdp"] <- "FDP (de)"
labs[labs=="social.DE.grun"] <- "Greens (de)"
labs[labs=="social.DE.linke"] <- "Linke (de)"
labs[labs=="social.DE.afg"] <- "AfD (de)"
#
labelsL <- paste(labs, "   ")
labelsR <- paste("   ", labs)
#
labelsL[11:21] <- ""
labelsR[1:10] <- ""
#
colors <- rep(viridis(21, begin=0.1, end=0.9, direction=1), 2)
dat2$colors <- rep(viridis(21, begin=0.1, end=0.9, direction=1), 2)
#
# Figure A1b. Results from alternate specifications of the joint Bayesian Aldrich-McKelvey model.
#
pdf("figA1b.pdf", height=9, width=8)
#
#
ggplot(dat2) +
  geom_line(aes(x = as.factor(variable), y = value, group = id, color=colors), size = 2) +
  geom_point(aes(x = as.factor(variable), y = value, color=colors), size = 5) +
  theme_minimal(base_size = 18) +
  scale_color_manual(values=c(dat2$colors[1:21])) +
  #scale_color_manual(values=rep(viridis(24, begin=0.1, end=0.9, direction=1), 4)) +
  #scale_color_viridis(discrete=TRUE, option="C", end=0.9, direction=1)   +
  #scale_color_distiller(palette = "Spectral") +
xlab("") +
geom_text(data = subset(dat2, variable=="No vignettes"),
		aes(x = as.factor(variable), y = value, color = colors, label = labelsL),
        size = 5, hjust = 1, vjust=c(0.0,  0.5,  0.5,  0.5,  0.5, -0.5, 1.5, -0.5,  0.0,  -0.75,  0.5, -0.5,  0.5,  0.5,  0.5,  0.0,  0.0,  0.5,  0.5,  0.0,  0.0)) +
geom_text(data = subset(dat2, variable=="Only DR"),
        aes(x = as.factor(variable), y = value, color = colors, label = labelsR),
        size = 5, hjust = 0, vjust=sample(c(-0.5,0,0.5),21,replace=TRUE)) +
  theme(legend.position = "none",
#      panel.grid.major.y = element_blank(),
      panel.grid.minor.y = element_blank(),
#      panel.grid.major.x = element_blank(),
      axis.ticks.y = element_blank(),
      axis.ticks.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.y = element_blank(),
      plot.title = element_text(size = 18, hjust=0.5)) +
      ggtitle("Social left-right estimates")
#
dev.off()
#
sink()
#
